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A new technique Is described for solving supersonic fluid dynamics problems containing multiple regions of 
continuous flow, each bounded by a permeable or Impermeable surface. Region boundaries are, in general, 
arbitrarily shaped and time dependent. Discretization of such a region for solution by conventional finite dif- 
ference procedures is accomplished using an elliptic solver which alleviates the dependence on a particular base 
coordinate system. Multiple regions are coupled together through the boundary conditions. The technique has 
been applied to a variety of problems including a shock diffraction problem and supersonic flow over a pointed 
ogive. 


Nomenclature 


G;, Cg 

E 


P 

f 

G 

G 

S 

I 

J 

Ms 

P 


P 

Q 

Q 

Q 

f 

S 

5/ 

5/ 


= coefficients used in matrix S 
= error norm in integrated Jacobian 
= total energy per unit volume 
= {-flux vector 
=jc-flux vector 
= grid operator 
= r;-flux vector 
=y-flux vector 
= inverse Jacobian 
= Jacobian 

= shock Mach number 

-{-forcing function for grid and grid speed 
equations 
= static pressure 

= 17 -forcing function for grid and grid speed 
equations 

= scaled solution vector, dependent variable 
vector 

== solution vector, dependent variable vector 
= right-hand side of linear grid speed equation 
= coefficient operator of linear grid speed 
equation 

= positive shift operator indicating a forward 
shift ony-index (i.e., Sf[<t>j^A^^j+i.k) 

= negative shift operator indicating a backward 
shift on y-index (i.e., Sj{4>j^A 
= positive shift operator indicating a forward 
shift on /:-index (i.e., ) 

= negative shift operator indicating a backward 
shift on A:-index (i.e., SjA^jk] = <t>ik-j) 


t 

u,v 

x,y 

z 

a,0 

y 

y 

A 

5 


p 

T 

V 

00 


Subscripts 

g 


fk 

t,X,y,ri,^,T 


iv 

Superscripts 

a,b 

n 


= time 

-x,y velocity components, respectively 
= Cartesian coordinates 
= grid speed vector 

= coefficients in grid and grid speed operators 
= T,{,Tj direction finite-difference operators, 
respectively 

= coefficient in grid and grid speed operators 
= ratio of specific heats 

= forward difference operator or denotes in- 
cremental value 
= central difference operator 
= upper bound for error norm in integrated 
Jacobian 

= general curvilinear coordinates 
= density 

= general curvilinear marching coordinate 
= backward difference operator 
= infinity 


= indicates association with actual finite- 
difference grid 

= indices in {,77 directions, respectively 
= partial differentiation with respect to this 
independent variable 

= second partial derivative with respect to this 
independent variable 
= cross partial derivative 


= used to identify two numerical represen- 
tations of the same metric quantity 
= marching direction index 
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Introduction 

T he electronic revolution, which undoubtedly is still in its 
infancy, has given man an incredibly powerful tool with 
which to solve many of the problems faced in today’s highly 
technological society. This paper deals with the application of 
this tool, the high-speed digital computer, to problem solving 
in fluid dynamics. 

Since the advent of the high-speed digital computer, an 
extensive effort has been made toward obtaining solutions to 
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fluid dynamics problems for which there exists a concise 
mathematical treatment resulting in a system of equations and 
boundary conditions describing an approximation to various 
physical processes. Many different mathematical descriptions 
exist depending upon the nature of the assumptions inherent 
to the derivation. The present work concentrates on fluid 
dynamics problems in which the flow is assumed to have a 
compressible, nonviscous nature. Such flow problems are 
adequately described by the Euler equations. An additional 
simplification (not restriction) which is made in the present 
work is that of two-dimensional flow. 

There are many computer codes in existence for solving the 
two-dimensional Euler equations. Unfortunately, each of 
them was written with a particular class of problems in mind. 
This places restrictions on the range of applicability of any 
computer code. Such restrictions may be categorized as those 
pertaining to the physics of the flow itself and those per- 
taining to the problem geometry description. These two 
categories are not necessarily independent. The former 
category includes such things as the development of flow 
singularities, discontinuities, or steep gradients for which 
there is inadequate numerical treatment. The latter category 
involves the manner in which the flowfield grid point 
distribution is defined relative to a base coordinate system. 
For example, a typical restriction of this type arises when a 
shock boundary is defined in Cartesian variables as 
Xg =x^ (yj) and possibly due to some interaction process, the 
shock slope, dx^/dy, becomes unbounded at some point. This 
difficulty can certainly be remedied in a given situation by a 
coordinate rotation or some other trick, but in the general 
case, any such trick, which may well require extensive code 
modification, will undoubtedly have similar limitations of its 
own. It is the purpose of this research effort to remove some 
of these restrictions by developing a generalized two- 
dimensional Euler equation solver using a modular approach 
and a very general treatment of the module geometry, thus 
providing one computer code capable of solving a wide 
variety of two-dimensional fluid dynamics problems. 

Modular Approach 

Although the present approach and the resulting computer 
code are in no way limited to supersonic flow problems, such 
problems do provide a more extensive test of the general 
concepts developed. As a result, the discussion presented 
herein will tend to emphasize the application of the present 
approach to supersonic flow. 

Many steady and unsteady supersonic flow problems 
contain multiple regions of continuous flow, each of which is 
either bounded by a surface across which the flow is 
discontinuous such as shock waves or slip surfaces or 
bounded by an impermeable surface such as the body. The 
present effort is to develop a cornputational solver which is 
designed to compute the solution to the Euler equations in an 
arbitrary region or module. The complete flow problem may 
consist of many of these modules coupled together through 
the appropriate application of boundary condition 
procedures. For example, the single Mach reflection planar 
shock diffraction problem shown in Fig. 1 may be described 
with two modules. These modules share a slip surface as a 
common boundary. 

Since the present approach requires that the flow module 
boundaries also be computational boundaries a generalized 
mapping of the independent variables must be performed (see 
Fig. 2). It is clear that in the general case these module 
boundaries are time-varying in nature. Such boundaries must 
therefore be capable of assuming virtually any shape dictated 
by the governing equations. Consequently, it is especially 
important that there is no built-in dependence of the validity 
of the module geometry description (grid point distribution 
and movement) on the particular base coordinate system 
chosen as a reference frame. The unique manner in which the 
present technique avoids such dependence is described in a 
later section. 



Physical Plane 
ShKk -Diffraction Problem 



Computational Plane 
Two -Solver -Net 

Fig. 1 Conceptual network of solvers. 


A literature search reveals that while extensive information 
exists on the idea of patching together of solutions, for 
example in boundary layer-inviscid interactions, etc., very 
little information is available on the present modular type of 
approach. Ludloff and Friedman^ solve the Mach reflection 
planar shock diffraction problem with what is apparently two 
modules although they do not specifically indicate such an 
approach. 

In contrast to the lack of information available on the 
modular type approach, considerable information is available 
on the subject of generalized geometry. Various types of 
automatic grid generation procedures have been developed ^ 
some of which have been applied to domains with moving 
boundaries. The present approach uses the automatic 
grid generation procedure of Thompson et al.^ and extends it 
in a unique way to allow for domains with time varying as 
well as arbitrarily shaped boundaries. 


Governing Equations 

The two-dimensional unsteady Euler equations are written 
in conservation-law form in Cartesian coordinates as 


dt dx dy 


( 1 ) 


where are given by 





" pu - 


pv ^ 


pu 

./= 

p^pu^ 


puv 

Q = 

pv 

puv 


p-^pv^ 





^{p^e)u^ 


^ {p + e)v^ 


with {u,u) representing the Cartesian {x,y) velocity com- 
ponents, p the density, p the static pressure, and e the total 
energy per unit volume, e is related to p»p,u,v through the 
equation 


e=-^ + ^{u^ + v^) 

'V — / 2 

where ^ is the ratio of thermal capacities of the fluid. 

Cartesian coordinates are used as the base coordinate 
system but in order to map bodies and other surfaces onto 
constant coordinate lines, the following coordinate mapping 
is introduced (see Fig. 2): 

r = = (2) 
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Computational Plane 

Fig. 2 Transformation of physical plane to computational plane for 
single solver. 


The governing equations, Eq. (I), are transformed by this 
mapping according to Viviand*® into a new strong con- 
servation-law form written as 


dQ dF dG 
— - 1 - — -0 
dT dr} 

where 

Q=Iq 

f'=H^,q + ^J+^yg) 

G = HVlQ+Vxf+Vyg) 

and 


( 3 ) 


( 4 ) 


where J is the Jacobian of the mapping and / is the Jacobian 
of the inverse mapping. The metrics of the mapping are 
related to the metrics of the inverse mapping by the equations 

Jr}^^-y^, IVy=x^ (5) 

Using Eq. (5), Eq. (4) is rewritten as 

{y,x^ -x,y^)i)-hyj-x^g 
( 5 = (x,y^ -y,x^ )Q-yJ-\‘X^g ( 6 ) 


Geometry 

There appears to be no universal procedure in the literature 
for treating the metric and Jacobian terms appearing in Eq. 
(6). This section, therefore, attempts to identify some logical 
ground rules which are followed In the treatment of these 
terms. The concepts involved are more completely discussed 
by Steger. * The rules stem from accuracy considerations and 
are based on the intuitive suggestion that for a scheme to be 
considered acceptable, the flowfield code which it supports 
must be capable of exact reproduction of a uniform flow. 
That is, with the boundary values held fixed at some uniform 
flow conditions and the initial flowfield set to these same 
conditions, the finite difference algorithm should exactly 
reproduce this same flowfield for all time. This is actually a 
statement of the independence that should exist between the 
physical flow and the grid distribution for (he case of uniform 
flow. It provides a simple test for an existing numerical 
algorithm and a criterion which ties down many of the 
questions arising in code development work regarding the 
manner in which the metrics of Eq. (6) should be computed. 

Consider the expanded form of Eq. (3); 

dlq d 

-x^y,)Q+yJ-x^g] 

d 

+ ^ -yrX( ) q -yj+x^ g\=o (7) 

Replace d/dr, d/d^, dfdr} by some finite difference operators, 
say 


dr 




dr} 




7 


resulting in 

r. [/^l + [(y,jr‘ -x,y \ ) ^] + T, [ {x,y\ -y,x \ ) q] 

+ + r, [ ~y\J\ + Fj [ - jf J#] + r, [;rf = 0 (8) 

where the superscripts a and b appearing on the metric terms 
arc used for later reference to identify two different numerical 
representations of the same quantity. Now for a uniform 
flow, since f=f{q), g=g(q), and it is required that 
q = constant throughout the flowfield for all time, then 

r,[<7] = iq ] = r, [9] = Fj [/I = r, t/1 = f^ [g] = f, [g] = 0 


and the following conditions result; 


r{l)'?J-r,ty»]=0 (10) 

Fj[jrJ]-F,[jff] = 0 (11) 

Note that the differential analogs of these equations are 
simply identities for a well behaved mapping. Clearly, Eqs. 
(10) and (11) are satisfied if the metrics are differenced with 
the same operators as those used in the finite difference 
scheme. Equation (9), however, is a numerical representation 
of an identity coined by Thomas and Lombard’* as the 
geometric conservation law (GCL). It says that the GCL 
equation 


Conservation-law form of the governing equations is 
necessary according to Lax’*^ to ensure that the jump con- 
ditions existing across weak solutions are automatically 
satisfied. This form of the governing equations thus adds a 
shock capturing capability to the resulting computer code. 


dl d d 

j^ + -^lyrX,-x,y^]+-lx,y(-y,x^] = 0 ( 12 ) 

must be differenced in an identical manner to the flow 
equations, Eq. (7), This result is nontrivial only in the case of 
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a time varying grid. Thomas and Lombard'* reached the 
same conclusion with consistency arguments and analogy with 
finite volume methods. It must be pointed out, however, that 
conditions (9-1 1) insofar as the present effort is concerned are 
strictly a result of choosing the strong conservation law 
representation of the governing equations. If the weak 
conservation-law equations had been chosen, a different set 
of conditions would result, this time for the differencing of 
the metric gradients appearing in the source terms. The use of 
the nonconservative form of the equations results in no 
special geometry differencing requirements whatsoever. 

Thus far two sets of metric values have been identified. 
They are referred to as a-metrics and !?-metrics based on their 
superscript. Conditions (10) and (11) dictate the manner in 
which the !?-metrics are to be computed. However, the 
calculation of the a-metrics is still a free choice. Although the 
two sets represent the same physical quantities, they need not 
be numerically equivalent. In fact, it is shown in a subsequent 
section that the manner in which the a-metrics are computed is 
dictated by the accuracy with which the integrated Jacobian 
value, /, resulting from Eq. (9), represents the actual Jacobian 
of the mapping. 

The calculation of the metrics requires knowledge of the 
coordinates (x,y) of each grid point. Determination of these 
coordinates and the speed with which the points move {x^,y^) 
is the subject of the next section. 

Grid and Grid Speed Operators 

In order to determine the metric quantities, the coordinates 
{x,y) of each grid point must be known. Also, grid point 
speeds {x^,y ^ ) are required to advance both the flow solution 
and the Jacobian [Eqs. (7) and (12)) in time. Due to the time 
varying nature of the grid boundaries, the location and speed 
of each interior grid point are necessarily dependent upon the 
location and speed of the boundary points. Two methods for 
obtaining such dependence are now described. This elliptic- 
type dependence may be obtained for the grid by following the 
approach of Thompson et al.^ Given the boundary point 
coordinates {x,y), the interior grid point coordinates are 
required to satisfy the nonlinear elliptic coupled partial 
differential equations 


G[x] = 0, G[y]=0 

with specified boundary values, where 

G — ct — 2B 

(13) 



(14) 

where 



(15) 


(16) 

II 

(17) 


The forcing functions P(r,{,7;) and G(r,{,?y) are used to 
concentrate the grid lines where they are most needed. For the 
present study these functions are set to zero. 

The requirement that the coordinates (x,y) satisfy Eq. (13) 
plus boundary conditions allows determination of the grid. 
However, the grid speed values are still unknown on the 
interior points. The boundary point speeds are assumed 
known since these values typically represent the speed of 
shock points, etc., which are determined from the flow 
solution. Interior values for (x^,y^) could, of course, be 
obtained with backward finite differences but this requires 


extra information at the initial data surface and may be in- 
consistent with the scheme used to advance the boundary 
point locations in time. In addition, the solution to Eq. (13) 
must be iterative due to the nonlinearity of the operator, C. 
Another approach, and the one used in the present study, is to 
differentiate Eqs. (13) with respect to t which, for a one-to- 
one mapping, yields the equation 

S[Z]=r (18) 

where 

Z={x^,y^) 

r= -P ) (19) 


G + C,(t,{. ij) — 
Ok 

Ok 

+ Q(t.{.ii) ^ 

dri 

+ C<(t,{,j/) ^ 
On 

Ok 

PC 

On 

-^Cg(Tgi,7}) — - 
dri 

j 


and 

c, (T.k.r,) =2(x^Xf + /(Pjfj + Cat, )y^ ) 

Cj (T.i.ri) -I(PXf + Qx^ )y^ ) 

Cj (T.i.ri) -I(Px^ + Qjf, ) jr, ) 

(r.t.ij) =2{x^^y^ - W£ + Qx^)x ^ ) 

Cs (r.lr,) =2(y^Xf -y^^x^ +HPy( + Qv, )>-,) 

Ci(T.tv) =2(yux^ -liPy^ + Qy, )y( ) 
-y(^,-f(Py(+Qy,)x^) 

Cg{T.^.ti) =2(>{{3', -yi^yf +HPy( + Qy, )jfj ) (20) 

Once the coordinates (x-,>') are known at each grid point, 
the metric quantities and their derivatives may be obtained 
with finite differences. The result is that the system of partial 
differential equations represented by Eq. (18) is linear in the 
dependent variables (x^,y^) with known variable coefficients 
and a direct method of solution to its finite difference 
representation may be employed to determine (x^,y^) at all 
interior points when given the boundary point speeds. In 
addition, the grid-point locations may be determined from a 
simple time integration of these computed speeds rather than 
by solving the nonlinear system, Eq. (13). This point is 
discussed in the next section. 

Coupling of Geometry Treatment and 
Finite Difference Scheme 

Several points are considered in this section involving the 
accuracy of the procedures developed. This accuracy is in- 
timately connected to the coupling which exists between the 
finite difference scheme chosen to integrate the governing 
equations and the treatment of the geometry. MacCormack’s 
standard unsplit predictor-corrector scheme'’ is used to 
integrate the flow equations, Eq. (7), and the GCL equation, 
Eq. (12), in time. Application of this scheme yields the 
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corresponding finite-difference equations: 

Predictors 

(Iq) = (Iq) " - ^ A, I 0-,at; -x,y^)q+y^J-x^4] " 
Ar 


where fi denotes the central difference operator and S denotes 
the shift operator, then 

e^O{Ar^, Ar^Af, At^Ai;) 

This is a higher order than the accuracy of the numerical 
determination of from the formula 


„ Ar 


jn-¥l — 


A{ 




Ar 

— A* lx,yl-y,x’(]’‘ 


( 22 ) 


An additional result of obtaining (x.y)"*' from Eqs. (25) 
and (26) is that ^ 

IC"+'[Ar'’+^]l=e„ | = f^ 


Correctors 

(/<?)-' = ^ [(/<?) ^ + (/^) " - g V, [ (y,x‘ -x,y ^, ) q 

+yif-x^g]^'- p V* 

Aij 

X [ (x,y^ -yrXf)4-yf/+xlg] (23) 

/-+/ = ^ +/, _ 0 _^^y, j 


where f,j~0(AT^). That is {x,y) at the new time are very 
accurate representations of the values obtained by solving Eq. 
(13) at /I + 1. This fact provides some assurance that the grid 
will remain nicely structured as it moves in time. 

Algorithm 

The algorithm for applying the preceding procedures to an 
arbitrary module is described in the following 15 steps. A 
priori knowledge of the initial boundary point locations and 
their speeds, and the initial flow solution, is assumed. 

1) Given at all boundary points, compute 

at all interior points initially by solving the coupled equations 
CnA:'*]=0,C«l>'']=0. 

2) Compute the ^-metrics from Eq. (27). 

3) Compute the ^-metrics from 


- klx,y\-y^\]'<^i]^ (24) 


where A, V are the usual forward and backward finite- 
difference operators. 

The value of 7"+^ obtained from Eq. (24) must accurately 
reflect the true grid structure at « + 1. The grid must therefore 
be advanced in time in such a way that the Jacobian, 7J+^ as 
computed with finite differences from the new grief point 
locations, is an accurate representation of 7'’+^. 

That is 

£=I7;-^^-7«+/|:S€ 

where e is a small number. Two methods exist for controlling 
the order of c. The manner in which the new grid is obtained 
from a time integration controls the value of 7J+^ and the 
manner in which the j-metrics in Eqs. (21-24) are treated 
controls the value of 7"+^ It is shown by Hindman “ that if 
the grid is advanced in time by the Euler predictor-modified 
Euler corrector scheme. 


x''^’=x"+Atx^ 



y"*' =y +Ary” 


(25) 

X”*' = ‘/2 (x" +JC^ 

+ Atx?+') 



+ Ary'} * ') 

(26) 

and the j-metrics are differenced as (with A{ = At? = 1) 



yf=s,s^ty] 



yf=SjSf[y] 

(27) 


yf*'=6,S,^y^n 



yr=5^s;iy^n 

(28) 


yf xf =AJjc«] 

yf=Ajly^], Arf"=A,U«] 

4) Given (jr;,>;) at all boundary points, compute (x^,y^) 

at all interior points by solving Eq. (1 8). ^ ^ 

5) Compute the Jacobian 7" from the equation, 

r=p^=6j[xn^ lyn-hlxn^jly''] 

6) Apply Eqs. (21) and (22) to yield (7^) ^ and 7^^ at all 
interior grid points. (Special difference equations are used at 
boundary points.) 

7) Apply Eq. (^ to yield^"+^Z^0 at all grid points. 

8) Compute q"+'=(/q) »+'//"+' and _apply boundary 

conditions to obtain the boundary speed (x”*', y”+> ) (i . e . , 

shock speed, etc.). 

9) Compute a-metrics from Eq. (28). 

10) Compute b-metrics from 

yf*' = ^jly"*']. xf*‘ = Vj[x^'] 

U) Gi^n (x^'*^', y^*') at all boundary points, compute 
(x^*. y”*') at all interior points by solving Eq. (18). 

12) Apply Eqs. (23) and (24) to yield (Iq) ’•+>, I"*' at all 
interior grid points. (Special difference equations are used at 
boundary points.) 

13) Apply Eq. (26) to yield (x"*'.y+>) at all grid points. 

14) Compute q"*^ = (Ig) ”+> //"-n and apply boundary 

conditions to correct q"*’ and provide (x^*', y; + ') on the 
boundaries. ^ ^ 

15) Go to step 2. 


Initial and Boundary Conditions 
Initial Conditions 

The techniques developed in the present effort revolve 
around the concept of generality. The determination of an 
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Fig, 3 Grid for regular reflection planar shock diffraction, ramp 
angle = 60 deg, shock Mach number = 4.71. 



Fig. 4 Wall/ramp pressure distribution for regular reflection planar 
shock diffraction, ramp angle = 60 deg, shock Mach number = 4.71, 


initial flowfield for an arbitrary problem does not lend itself 
to a general treatment. As a result, the initial solution for a 
given flow configuration is treated as an independent problem 
and is not part of the existing computer code. 

Boundary Conditions 

The algorithm described in the previous section applies to 
an arbitrary module. Steps 8 and 14 require the application of 
boundary condition procedures. The space limitation prevents 
a discussion of these procedures here but such a discussion 
may be found in Ref. 21 . 

Numerical Results 

This section consists of two parts. The first deals with 
problems which were solved with one module. These 
problems provide an excellent test of the new geometry 
procedures developed in this study. The second part deals with 
a simple preliminary test of the multiple module capability for 
a problem with two modules and one interface boundary. 

Single Module 

Numerical results for two problems requiring a single 
module are presented. The results for a regular reflection 
planar shock diffraction problem with a shock Mach number 
of 4.71 and a ramp angle of 60 deg are shown in Figs. 3-5 and 
compared with the results of Kutler and Shankar pjgg 4 



Fig. 5 Wall/ramp density distribution for regular reflection planar 
shock diffraction, ramp angle = 60 deg, shock Mach number = 4.71 . 



Fig. 6 Grid for single module ogive with geometric singularity on 
shock, Mach number = 2.0, thickness = 0.1 chord. 


and 5. Both the pressure and density distributions agree well 
with the Kutler and Shankar solution. It is interesting to note 
the behavior of the pressure profile for this problem as it 
approaches the stagnation point. The profile tends toward a 
zero slope when approaching along the stagnation streamline 
but tends to spike as the stagnation point is approached along 
the ramp. This problem is unique in that the grid never 
reaches a steady-state solution. As 00, ^x/t and y ~*y/T 
due to the self-similarity that exists with respect to timk The 
meaning of a converged grid in this case is when tx^-x/rl 
and ly^ -y/rl are sufficiently small. Based on this definition, 
the converged grid is shown in Fig. 3. Note that a vortical 
singularity exists for this problem and is clearly shown in Fig. 
5. This singularity was captured by the present finite dif- 
ference algorithm. The computational domain for this 
problem represents a four-sided physical region. The 
supersonic inflow side of this region is an artificial boundary 
and its validity is due to the two-dimensional flow region 
which exists between the sonic line, shown in Fig. 3, and the 
shock ramp intersection point. If the numerical solution were 
required clear up to this intersection point, then the physical 
domain would have only three sides. The next example 
illustrates that such three-sided regions pose no serious threat 
to the present method. 

Consider a 10% thick pointed ogive body immersed in a 
Mach 2 freestream. The converged grid is shown in Fig. 6, 
thus illustrating a three-sided physical domain. This region is 
mapped into a computational four-sided region by in- 
troducing a “corner” point along the leading edge shock. The 
Jacobian of the inverse transformation vanishes at this ar- 
tificial corner thus creating a geometric singularity. This type 
of singularity is simply an additional corner point option in 
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Fig. 7 Surface pressure coefficient for single module ogive with 
geometric singularity on shock, Mach number =2.0, thickness = 0.1 
chord. 



Fig. 8 Grid for double module ogive with aft region geometric 
singularity, Mach number = 2.0, thickness = 0,1 chord. 


the existing computer code. This option is explained by 
Hindman.^® Solutions were also obtained with this singularity 
introduced on the ogive surface and on the supersonic outflow 
boundary instead of on the shock. No difficulties were en- 
countered with any of these cases. Figure 7 illustrates the 
ogive body surface pressure distribution along with the 
numerical results of Schiff,^^ The comparison is good. 
Special attention must now be directed to the grid structure in 
the vicinity of the geometric singularity. The grid is extremely 
nonorthogonal in this region. Some researchers have 
suggested in recent times that such nonorthogonality is un- 
desirable from the standpoint of causing numerical dif- 
ficulties or inaccuracies in the solution. This is certainly not 
true for the pointed ogive problem presented here. 

Two Modules 

The ogive body just discussed also serves as a test case of a 
double module problem with the trailing edge shock forming 
the interface boundary between the two modules. Figure 8 
illustrates the converged grid and both the leading edge and 
trailing edge shocks. In this case, the aft region is three-sided 
and a geometric singularity is introduced along the symmetry 
boundary. The surface pressure distribution along the ogive 
and the symmetry boundary is depicted in Fig. 9 and com- 
pared to the shock capturing results of Schiff.^^ Generally 
good agreement is observed where the comparison can be 
made. The present solution, however, exhibits a slightly 
higher estimate of the pressure just behind the trailing edge 
shock. This reflects the fact that the shock slope is slightly 



Nondimensional Distance Along Body and 
Aft Symmetry Boundary 


Fig. 9 Surface and aft symmetry boundary pressure coefficient for 
double module ogive with aft region geometric singularity, Mach 
number = 2.0, thickness = 0.1 chord. 


larger than it should be at this point. This same behavior is 
observed behind the leading edge shock and is due to some 
slight inconsistency in the scheme used to treat this type of 
corner point. The anomaly appearing near the center of the 
aft region symmetry boundary is due to the existence of the 
geometric singularity. This is contrary to the findings for the 
three-sided region in the single module ogive case where no 
anomaly appeared. A possible explanation of this follows. 
The grid for the single module case (Fig. 6) appears to be 
nearly symmetric in some sense with respect to the singular 
point. Therefore, the truncation error in the solution, which is 
weakly dependent upon the problem geometry, will be 
essentially the same on either side of this point. Thus the 
solution will not be forced by the geometric contribution to 
the truncation error to exhibit a nonsmooth behavior through 
the singularity. On the other hand, this symmetry does not 
exist for the double module case as seen in Fig. 8. This causes 
the truncation error on one side of the singular point to be 
different from that on the other side which givei rise to a weak 
perturbation on the solution at points near the singularity 
where the geometric variables are rapidly changing. 

Concluding Remarks 

A general method is presented for solving the unsteady two- 
dimensional Euler equations on multiple flow regions with 
arbitrarily-shaped and time-varying boundaries. The method 
is applicable to problems with moving boundaries provided the 
velocity of such movement can be determined or specified. 
This includes problems with moving pistons, structural 
deformations, accelerating bodies, moving or stationary 
discontinuity surfaces such as shocks and slip surfaces, etc. In 
the case of discontinuity surfaces, the scheme has the 
capability of capturing any discontinuities whose approximate 
shape and location is not known a priori provided the strength 
of such discontinuities is not excessive. 

The resulting computer code may be used to solve a wide 
variety of two-dimensional flow problems. The multiple flow 
region capability may be used to compute one flowfield with 
multiple regions or it may be used to simultaneously perform 
a parametric study of the solution for one flow problem. 
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